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AN APPROACH TO ASSESS DELAMINATION PROPAGATION SIMULATION 
CAPABILITIES IN COMMERCIAL FINITE ELEMENT CODES 

Ronald Krueger' 

ABSTRACT 

An approach for assessing the delamination propagation simulation capabilities in 
commercial finite element codes is presented and demonstrated. For this investigation, the 
Double Cantilever Beam (DCB) specimen and the Single Leg Bending (SLB) specimen 
were chosen for full three-dimensional finite element simulations. First, benchmark results 
were created for both specimens. Second, starting from an initially straight front, the 
delamination was allowed to propagate. The load-displacement relationship and the total 
strain energy obtained from the propagation analysis results and the benchmark results were 
compared and good agreements could be achieved by selecting the appropriate input 
parameters. Selecting the appropriate input parameters, however, was not straightforward 
and often required an iterative procedure. Qualitatively, the delamination front computed 
for the DCB specimen did not take the shape of a curved front as expected. However, the 
analysis of the SLB specimen yielded a curved front as was expected from the distribution 
of the energy release rate and the failure index across the width of the specimen. Overall, the 
results are encouraging but further assessment on a structural level is required. 


1. INTRODUCTION 

One of the most common failure modes for composite structures is delamination [1-4]. The 
remote loadings applied to composite components are typically resolved into interlaminar tension 
and shear stresses at discontinuities that create mixed-mode I, II and III delaminations. To 
characterize the onset and propagation of these delaminations, the use of fracture mechanics has 
become common practice over the past two decades [1, 5, 6], The total strain energy release rate, 
Gr, the mode I component due to interlaminar tension, Gi, the mode II component due to 
interlaminar sliding shear, Gn, and the mode III component, Gui, due to interlaminar scissoring 
shear, as shown in Figure 1, need to be calculated. In order to predict delamination onset or 
propagation for two-dimensional problems, these calculated G components are compared to 
interlaminar fracture toughness properties measured over a range from pure mode I loading to pure 
mode II loading [7-9]. A quasi static mixed-mode fracture criterion is detennined by plotting the 
interlaminar fracture toughness, G, , versus the mixed-mode ratio, Gii/Gt, detennined from data 
generated using pure mode I Double Cantilever Beam (DCB) (Gn/G-/=0), pure mode II End- 
Notched Flexure (ENF) (Gu/Gt=\), and Mixed-Mode Bending (MMB) tests of varying ratios, as 
shown in Figure 2a for T300/914C and Figure 2b for C12K/R6376 [10, 11]. A curve fit of these 
data is performed to detennine a mathematical relationship between G c and Gu/Gt. [12, 13], Failure 
is expected when, for a given mixed-mode ratio Gu/Gt, the calculated total energy release rate, Gt, 
exceeds the interlaminar fracture toughness, G c . An interaction criterion incorporating the scissoring 
shear (mode III), was recently proposed by Reeder [14]. The edge-cracked torsion test (ECT) to 
measure Gm c is being considered for standardization [15, 16]. 


1 National Institute of Aerospace (NIA), 100 Exploration Way, Hampton, VA 23666, resident at Durability, Damage 
Tolerance and Reliability Branch, NASA Langley Research Center, MS 188E, Hampton, VA, 23681. 
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The virtual crack closure technique (VCCT) is widely used for computing energy release 
rates based on results from continuum (2D) and solid (3D) finite element analyses and to supply the 
mode separation required when using the mixed-mode fracture criterion [17, 18]. The virtual crack 
closure technique has been used mainly by scientists in universities, research institutions and 
government laboratories and is usually implemented in their own specialized codes or used in post- 
processing routines in conjunction with general purpose finite element codes. An increased interest 
in using a fracture mechanics based approach to assess the damage tolerance of composite structures 
in the design phase and during certification has also renewed the interest in the virtual crack closure 
technique. The VCCT technique was recently implemented into the commercial finite element 
codes ABAQUS® 1 , NASTRAN" 2 and Marc™ 3 [19-21], The implementation into the commercial 
finite element code SAMCEF™ 4 [22] is a mix of VCCT and the Virtual Crack Extension Method 
suggested by Parks [23]. As new approaches for analyzing composite delamination are incorporated 
in finite element codes, the need for comparison and benchmarking becomes important. 

The objective of this study was to create an approach, independent of the analysis software 
used, which allows the assessment of delamination propagation simulation capabilities in 
commercial finite element codes. For this investigation, the Double Cantilever Beam (DCB) 
specimen with a unidirectional and a multi-directional layup and the Single Leg Bending (SLB) 
specimen with a multi-directional layup (as shown in Figure 3) were chosen for full three- 
dimensional finite element simulations. These specimen configurations were chosen, since they are 
simple and a number of combined experimental and numerical studies had been performed 
previously where the critical strain energy release rates were evaluated [24-27]. To avoid 
unnecessary complications, experimental anomalies such as fiber bridging were not addressed. 
Comparisons with test results will follow later in another report. First, benchmark results were 
created using models simulating specimens with different delamination lengths. For each 
delamination length modeled, the load and displacement at the load point were monitored. The 
mixed-mode strain energy release rate components were calculated along the delamination front 
across the width of the specimen. A failure index was calculated by correlating the results with the 
mixed-mode failure criterion of the graphite/epoxy material. It was assumed that the delamination 
propagated when the failure index reached unity. Thus, critical loads and critical displacements for 
delamination onset were calculated for each delamination length modeled. These critical 
load/displacement results were used as a benchmark. The computed total strain energy was also 
used as a benchmark. Second, starting from an initially straight front, the delamination was allowed 
to propagate based on the algorithms implemented into the commercial finite element software. The 
approach was demonstrated for the commercial finite element code ABAQUS 8 with focus on their 
implementation of the Virtual Crack Closure Technique (VCCT) [19]. VCCT control parameters 
were varied to study the effect on the computed load-displacement behavior during propagation. It 
was assumed that the computed load-displacement relationship should closely match the benchmark 
results established earlier. As a qualitative assessment, the shape of the computed delamination 
fronts was also compared to photographs of failed specimens. 


1 ABAQUS r is manufactured by Dassault Systemes Simulia Corp. (DSS), Providence, RI, USA 

2 NASTRAN®' is a registered trademark of NASA and manufactured by MSC. Software Corp., Santa Ana, CA, USA 

3 Marc™ and Mentat™ are manufactured by MSC. Software Corp., Santa Ana, CA, USA 

4 SAMCEF™ is manufactured by Samtech, Liege, Belgium 


2 



2. SPECIMEN DESCRIPTION 

For the current numerical investigation, the Double Cantilever Beam (DCB) and the Single 
Leg Bending (SLB) specimens, as shown in Figure 3, were chosen. The DCB specimen is used to 
detennine the mode I interlaminar fracture toughness, GY (Gn/Gf=0) [7]. The SLB specimen was 
introduced for the determination of fracture toughness as a function of mixed-mode I/II ratio [26, 
28], This test may be performed in a standard three-point-bending fixture such as that used for the 
ENF test. By varying the relative thickness of the delaminated regions (fi and ti), a different mixed- 
mode ratio may be achieved. This type of specimen was chosen to study mode separation. 
Previously, a number of combined experimental and numerical studies of these specimens had been 
performed and the critical strain energy release rates were evaluated [24-27]. 

In general, mode I, mode II and mixed-mode tests are performed on unidirectionally 
reinforced laminates, which means that delamination propagation occurs at a [0/0] interface and 
crack propagation is parallel to the fibers. For the current study, a DCB specimen made of 
T300/1076 graphite/epoxy with a unidirectional layup, [0] 2 4, was modeled. Although this 
unidirectional layup is desired for standard test methods to generate fracture toughness data, 
delamination propagation between layers of the same orientation will rarely occur in real structures. 
Previously, combined experimental and numerical studies on specimens with multi-directional 
layups were performed where the critical strain energy release rates of various interfaces were 
evaluated under mode I, mode II and mixed-mode conditions [25, 26]. Therefore, a DCB-specimen 
made of C12K/R6376 graphite/epoxy with a multi-directional layup was selected. The stacking 
sequence [±30/0/-30/0/30/0 4 /30/0/-30/0/-30/30/ T -30/30/0/30/0/ -30/0 4 /-30/0/30/0/±30] was 
designated D±30, where the arrow Q denotes the location of the delamination. Additionally, a SLB 
specimen with D±30 layup was also modeled. The material properties are given in Table I. 


3. METHODOLOGY 
3.1 Fracture Criteria 

Linear elastic fracture mechanics has proven useful for characterizing the onset and 
propagation of delaminations in composite laminates [5, 6]. When using fracture mechanics, the 
total strain energy release rate, Gt, is calculated along the delamination front. The term, GY, consists 
of three individual components, as shown in Figure 1. The first component, G/, arises due to 
interlaminar tension. The second component, Gu, arises due to interlaminar sliding shear (shear 
stresses parallel to the plane of delamination and perpendicular to the delamination front). The third 
component, Gm, arises due to interlaminar scissoring shear (shear stresses parallel to the plane of 
delamination and parallel to the delamination front). The calculated G/, Gn, and Gm components are 
then compared to interlaminar fracture toughness values in order to predict delamination onset and 
propagation. The interlaminar fracture toughness values are determined experimentally over a range 
of mode mixity from pure mode I loading to pure mode II loading [7-9]. 

A quasi static mixed-mode fracture criterion is determined by plotting the interlaminar 
fracture toughness, G c , versus the mixed-mode ratio, Gu/Gt. The fracture criteria is generated 
experimentally using pure Mode I ( G//G-/K)) Double Cantilever Beam (DCB) tests [7], pure Mode 
II (Gn/Gf= 1) End-Notched Flexure (ENF) tests [9], and Mixed Mode Bending (MMB) tests of 
varying ratios of G/ and Gu [8]. Typical examples are presented in Figure 2 for T300/914C and 
C12K/R6376 carbon epoxy materials. A 2D fracture criterion was suggested by Benzeggah and 
Kenane [13] using a simple mathematical relationship between G c and G///GY 


3 



77 


( 1 ) 


G c = G Ic + { G IIc~ G Ic)' 



In this expression, Gu and Gn c are the experimentally-detennined fracture toughness data for 
mode I and II as shown in Figure 2. The factor r/ was determined by a curve fit using the 
Levenberg-Marquardt algorithm in KaleidaGraph™ graphing and data analysis software [29]. 
Fracture initiation is expected when, for a given mixed-mode ratio Gn/Gp, the calculated total 
energy release rate, GY, exceeds the interlaminar fracture toughness, G c and therefore the failure 
index Gp/G c is equal or greater than unity 


G T 

~G C 


> 1 . 


( 2 ) 


For three-dimensional analysis, which yields results for the scissoring mode Gm, a modified 
definition is introduced where G s denotes the sum of the in-plane shearing components G n +G n i 
[30]. This modification becomes necessary if a mixed-mode failure criterion, which accounts for all 
three modes, is not available. For analyses where Gm= 0, this definition is equal to the commonly 
used definition of the mixed-mode ratio, Gjj/Gt mentioned above. To detennine failure along the 
delamination front, the critical energy release rate G c is calculated using equation (1) with Gu = Gs 
at each point along the delamination front. Subsequently, the failure index Gp/G c is detennined as 
above. The modified interaction criterion is an integral part of the VCCT for ABAQUS® analysis 
software and may be selected by the user [19]. 

Recently, Reeder [14] suggested an interaction criterion that is based on the 2D fracture 
criterion suggested by Benzeggah and Kenane [13] but incorporates the mode III scissoring shear 


G c = G Ic + ( G IIc - G Ic ) ' 


G I1 + G III 


+ { G IIIc- G IIc)' 


Gu + Gm 


G II + G III 


(3) 


which is also an integral part of the VCCT for ABAQUS 8 analysis software and may be selected by 
the user [19]. 

Although several specimens have been suggested for the measurement of the mode III 
interlaminar fracture toughness property [15, 31, 32] a standard does not yet exist. Currently, the 
edge-cracked torsion test (ECT) is being considered for standardization as a pure mode III test [15, 
16]. 


3.2 Virtual Crack Closure Technique (VCCT) 

3.2.1 Background 

A variety of methods are used in the literature to compute the strain energy release rate 
based on results obtained from finite element analysis. For delaminations in laminated composite 
materials where the failure is highly dependent on the mixed-mode ratio (as shown in Figure 2), the 
virtual crack closure technique (VCCT) [17, 18] has been most widely used for computing energy 
release rates. VCCT calculations using continuum (2D) and solid (3D) finite element analyses 
provide the mode separation required when using the mixed-mode fracture criterion. 
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The mode I, and mode II components of the strain energy release rate, G/, Gu are computed 
using VCCT as shown in Figure 4a for a 2D four-node element. The terms F ’ xi , F’ yi are the forces at 
the crack tip at nodal point i and u, , v\ and u' p , v[* are the displacements at the corresponding 
nodal points / and i behind the crack tip. Note that Gui is identical to zero in the 2D case. For 
geometrically nonlinear analysis where large deformations may occur, both forces and 
displacements obtained in the global coordinate system need to be transfonned into a local 
coordinate system (pc', v') which originates at the crack tip as shown in Figure 4a. The local crack tip 
system defines the tangential (pc', or mode II) and nonnal ( y ', or mode I) coordinate directions at the 
crack tip in the deformed configuration. The extension to 3D is straight forward as shown in 
Figure 4b and the total energy release rate Gt is calculated from the individual mode components as 
Gr =Gi +Gn +Gni- For the two-dimensional case shown in Figure 4a, Gui =0. 

3.2.2 VCCT for ABAQUS® 

Currently, VCCT for ABAQUS® is an add-on capability to ABAQUS “/Standard Versions 
6.5, 6.6 and 6.7 that provides a specific implementation of the virtual crack closure technique within 
ABAQUS 8 . The implementation of VCCT enables ABAQUS 8 to solve delamination and 
debonding problems in composite materials. The implementation is compatible with all the features 
in ABAQUS 8 such as large-scale nonlinear, models of composite structures including continuum 
shells, composite materials, cohesive elements, buckling, and contact. The plane of delamination in 
three-dimensional analyses is modeled using the existing ABAQUS “/Standard crack propagation 
capability based on the contact pair capability [19]. Additional element definitions are not required, 
and the underlying finite element mesh and model does not have to be modified [19]. 

Beyond simple calculations of the mixed-mode strain energy release rates along the 
delamination front, which was studied previously [27], the implementation also offers a crack 
propagation capability in ABAQUS 8 . It is implied that the energy release rate at the crack tip is 
calculated at the end of a converged increment. Once the energy release rate exceeds the critical 
strain energy release rate (including the user-specified mixed-mode criteria as shown in Figure 2), 
the node at the crack tip is released in the following increment, which allows the crack to propagate. 
To avoid sudden loss of stability when the crack tip is propagated, the force at the crack tip before 
advance is released gradually during succeeding increments in such a way that the force is brought 
to zero no later than the time at which the next node along the crack path begins to open [19]. 

In addition to the mixed-mode fracture criterion, VCCT for ABAQUS 8 requires additional 
input for the propagation analysis. If a user specified release tolerance is exceeded in an increment 
( G-G c )/G c > release tolerance, a cutback operation is perfonned which reduces the time 
increment. In the new smaller increment, the strain energy release rates are recalculated and 
compared to the user specified release tolerance. The cutback reduces the degree of overshoot and 
improves the accuracy of the local solution [19]. A release tolerance of 0.2 is suggested in the 
handbook [19]. 

To help overcome convergence issues during the propagation analysis, ABAQUS 8 
provides: 

• contact stabilization which is applied across only selected contact pairs and used 
to control the motion of two contact pairs while they approach each other in multi-body contact. 
The damping is applied when bonded contact pairs debond and move away from each other [19] 

• automatic or static stabilization which is applied to the motion of the entire model 
and is commonly used in models that exhibit statically unstable behavior such as buckling [19] 


5 



• viscous regularization which is applied only to nodes on contact pairs that have 
just debonded. The viscous regularization damping causes the tangent stiffness matrix of the 
softening material to be positive for sufficiently small time increments. Viscous regularization 
damping in VCCT for ABAQUS® is similar to the viscous regularization damping provided for 
cohesive elements and the concrete material model in ABAQUS®/Standard [19]. 

Setting the value of the input parameters correctly is often an iterative procedure, which 
will be discussed later. 


4. FINITE ELEMENT MODELING 

Typical three-dimensional finite element models of Double Cantilever Beam (DCB) and 
Single Leg Bending (SLB) specimens are shown in Figures 5 to 10. Along the length, all models 
were divided into different sections with different mesh refinement. A refined mesh of length 
d=5 mm with 20 elements was used for the DCB specimen as shown in the detail of Figure 5a. This 
section length had been selected in previous studies [24, 27] and was also used during the current 
investigation. Across the width, the model was divided into a center section and a refined edge 
section, j, to capture local edge effects and steep gradients. These sections appear as dark areas in 
the full view of the specimen as shown in Figure 5a. The specimen was modeled with solid brick 
elements C3D8I which had yielded excellent results in a previous study [27]. The DCB specimen 
with unidirectional layup, [0] 2 4, was modeled with six elements through the specimen thickness (2h) 
as shown in the detail of Figure 5 a. This model was used to calculate mode I energy release rates 
and create the benchmark results discussed later. For all the analyses perfonned, the nonlinear 
solution option in ABAQUS®/Standard was used. For propagation analyses using VCCT for 
ABAQUS 81 , the model with a unifonn mesh across the width, as shown in Figure 5b, was used to 
avoid potential problems at the transition between the coarse and very fine mesh near the edges of 
the specimen. 

For the analysis with VCCT for ABAQUS® 1 , the plane of delamination was modeled as a 
discrete discontinuity in the center of the specimen. To create the discrete discontinuity, each model 
was created from separate meshes for the upper and lower part of the specimens with identical nodal 
point coordinates in the plane of delamination [19]. Two surfaces (top and bottom surface) were 
created on the meshes as shown in Figure 5. Additionally, a node set was created to identify the 
intact (bonded nodes) region. Two coarser meshes with a reduced number of elements in width and 
length directions were also generated as shown in Figures 6a and b. 

Three models of the DCB specimen were generated with continuum shell elements SC8R as 
shown in Figures 7a to c. The continuum shell elements in ABAQUS ® are used to model an entire 
three-dimensional body, unlike conventional shells which discretize a reference surface. The SC8R 
elements have displacement degrees of freedom only, use linear interpolation, and allow finite 
membrane defonnation and large rotations and, therefore, are suitable for nonlinear geometric 
analysis. The continuum shell elements are based on first-order layer-wise composite theory and 
include the effects of transverse shear defonnation and thickness change [33]. In the x-y plane, the 
models have the same fidelity as the models made of solid brick elements C3D8I shown in Figures 
5b, 6a and 6b. In the z-direction, only one element was used to model the thickness of the specimen. 
These less refined models were used to study the effect on perfonnance (CPU time), computed 
load/displacement behavior and delamination front shape in comparison with the more refined 
model discussed above. 
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For the DCB specimen with multi-directional layup, D±30, a model with a unifonn mesh 
across the width was used as shown in Figure 8. The DCB specimen was modeled with solid brick 
elements C3D8I which had yielded excellent results in a previous study [27]. Two plies on each 
side of the delamination were modeled individually using one element for each ply as shown in the 
detail of Figure 8. Since the delamination occurs at an interface between materials with dissimilar 
properties, care must be exercised in interpreting the values for Gi and Gn obtained using the virtual 
crack closure technique. For interfacial delaminations between two differing orthotropic solids, the 
observed oscillatory singularity at the crack tip becomes an issue for small element lengths [34, 35]. 
Hence, a value of crack tip element length, A a, was chosen (approximately three ply thicknesses) in 
the range over which the strain energy release rate components exhibit a reduced sensitivity to the 
value of A a. The adjacent four plies were modeled by one element with material properties smeared 
using the rule of mixtures [36, 37], Smearing appeared suitable to reduce the model size, however, it 
did not calculate the full A-B-D stiffness matrix contributions of the plies. The adjacent element 
extended over the four 0° plies. The six outennost plies were modeled by one element with smeared 
material properties. 

For the SLB specimen with multi-directional layup, D±30, a model with a uniform mesh 
across the width was used as shown in Figure 9. The SLB specimen was modeled with solid brick 
elements C3D8I which had yielded excellent results in a previous study [27]. For modeling 
convenience, the upper and lower anns were modeled similar to the model of the DCB specimen. 
To model the test correctly, only the upper ann was supported in the analysis as shown in Figure 9. 
An additional mesh with a longer refined center section was generated as shown in Figure 10. The 
refined model was used to study the effect on computed load/displacement behavior and 
delamination front shape in comparison with the model discussed above. 


5. ANALYSIS 

First, models simulating specimens with different delamination lengths were analyzed. For 

each delamination length modeled, the load and displacement at the load point were monitored. The 
mixed-mode strain energy release rate components were calculated along the delamination front 
across the width of the specimen. A failure index was calculated by correlating the results with the 
mixed-mode failure criterion of the graphite/epoxy material. It was assumed that the delamination 
propagated when the failure index reached a value of unity. Thus, critical loads and critical 
displacements for delamination onset were calculated for each delamination length modeled. These 
critical load/displacement results were used as a benchmark. Second, starting from an initially 
straight front, the delamination was allowed to propagate based on the algorithm implemented into 
VCCT for ABAQUS®. Input parameters were varied to study the effect on the computed load- 
displacement behavior during propagation. It was assumed that the computed load-displacement 
relationship should closely match the benchmark results established earlier. 

The total strain energy in the model was calculated from the computed load/displacement 
behavior. The results were compared with the values computed internally by ABAQUS ® . The total 
strain energy was also compared to the damping energies associated with the different stabilization 
techniques in ABAQUS 5 * * 8 . Input parameters were varied to study the ratio between the damping 
energies and the total strain energy. It was assumed that input parameters which produced results 
with the smallest damping energies corresponded to results which also matched the benchmark 
results best. 
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As a qualitative assessment, the shape of the computed delamination fronts were also 
compared to photographs of failed specimens. 


5,1 Creating a Benchmark Solution for a DCB specimen with unidirectional layup 

The computed mode I strain energy release rate values were plotted versus the normalized 
width, y/B, of the specimen as shown in Figure 11. The results were obtained from models shown in 
Figure 5 a for seven different delamination lengths a. An opening displacement 6/2=1. Omm was 
applied to each ann of the model. Qualitatively, the mode I strain energy release rate is fairly 
constant in the center part of the specimen and drops progressively towards the edges. This 
distribution will cause the initial straight front to grow into a curved front as explained in detail in 
the literature [38-41], As expected, the mode II and mode III strain energy release rates were 
computed to be nearly zero and hence are not shown. Computed mode I strain energy release rates 
decreased with increasing delamination length a. 

The failure index Gj/G c was computed based on a mode I fracture toughness G/ r = 1 70.3 J/nT 
for T300/914C (see Figure 2a). The failure index was plotted versus the nonnalized width, y/B, of 
the specimen as shown in Figure 12. For all delamination lengths modeled, except for <7=40 mm, the 
failure index in the center of the specimen (y/B= 0) is above unity (Gr/G,>\). 

For all delamination lengths modeled, the reaction loads P at the location of the applied 
displacement were calculated and plotted versus the applied opening displacement 6/2 as shown in 
Figure 13. The critical load, P cr u, when the failure index in the center of the specimen (y/B= 0) 
reaches unity {Gp/G c ~ 1), can be calculated based on the relationship between load P and the energy 
release rate G [42], 


P 2 dC P 
2 dA 


( 4 ) 


In equation (4), Cp is the compliance of the specimen and dA is the increase in surface area 
corresponding to an incremental increase in load or displacement at fracture. The critical load P crit 
and critical displacement 6 cr iJ2 were calculated for each delamination length modeled 



and the results were included in the load/displacement plots as shown in Figure 14 (solid red 
circles). The results indicate that, with increasing delamination length, less load is required to extend 
the delamination. This means that the DCB specimen exhibits unstable delamination propagation 
under load control. Therefore, prescribed opening displacements 6/2 were applied in the analysis 
instead of nodal point loads P to avoid problems with numerical stability of the analysis. It was 
assumed that the critical load/displacement results can be used as a benchmark. For the 
delamination propagation, therefore, the load/displacement results obtained from the model of a 
DCB specimen with an initially straight delamination of <7=30 mm length should correspond to the 
critical load/displacement path (solid red line) in Figure 14. 
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5.2 Delamination Propagation in a DCB Specimen with Unidirectional Layup Using VCCT 
for ABAQUS® 

5.2.1 Computed load/displacement behavior for different input parameters 

The propagation analysis was perfonned in two steps using the model shown in Figure 5b 
for a delamination length 30 mm. In the first step, a prescribed displacement {812= 0.74 mm) was 
applied in two increments which equaled nearly the critical tip opening {8 cn J2= 0.75 mm) 
detennined in the analysis above for a delamination length of <7=30 mm. Dividing the first step into 
just two increments was possible, since the load-displacement behavior of the specimen up to failure 
was linear as shown in Figure 14. In the second step, the total prescribed displacement was 
increased {812= 2.8 mm). Automatic incrementation was used with a small increment size at the 
beginning (10 4 of the total increment) and a very small minimum allowed increment (10" 18 of the 
total increment) to reduce the risk of numerical instability and early tennination of the analysis. The 
analysis was limited to 1000 increments. Initially, analyses were perfonned without stabilization or 
viscous regularization. Release tolerance values between 0.2 and 0.6 were used. Using these 
parameters, the analysis tenninated early prior to advancing the delamination. 

In Figures 15 to 20, the computed resultant force (load P) at the tip of the DCB specimen is 
plotted versus the applied crack tip opening {8/2) for different input parameters which are listed in 
Table II. For the results shown, the analysis tenninated when the 1000 increment limit set for the 
analysis was reached. Several analyses tenninated early because of convergence problems. To 
overcome the convergence problems, the methods implemented in ABAQUS® were used 
individually to study the effects. For the results plotted in Figure 15, global stabilization was added 
to the analysis. For a stabilization factor of 2xl0" 5 , the stiffness changed to almost infinity once the 
critical load was reached causing the load to increase sharply (plotted in blue). The load increased 
until a point was reached where the delamination propagation started and the load gradually 
decreased following a saw tooth curve with local rising and declining segments. The gradual load 
decrease followed the same trend as the benchmark curve (in grey) but is shifted toward higher 
loads. For a stabilization factor of 2x1 O’ 6 (in green), the same saw tooth pattern was observed but 
the average curve was in good agreement with the benchmark result. For a stabilization factor of 
2x1 0' 7 (in red), the average was lower than before but was in good agreement with the benchmark 
result until tennination after 550 increments due to convergence problems. The results obtained for 
a stabilization factor of 2x1 0 s (in black for a release tolerance of 0.2) were on top of the previous 
result. The rate of convergence appeared to be slower since only 8/2= 1.14 mm was applied for 
1000 increments compared to 8/2= 1.24 mm for a stabilization factor of 2xl0" 6 and the same release 
tolerance (0.2). Changing the release tolerance also appeared to influence the convergence as shown 
in Table II. For a release tolerance of 0.02, the analysis tenninated after 1000 increments for 8/2= 
1.04 mm. For a release tolerance of 0.002, the analysis tenninated due to convergence problems 
after 45 1 increments. Changing the release tolerance, however, appeared to have no effect on the 
overall load/displacement behavior or the magnitude of the saw tooth pattern. 

An adaptive automatic stabilization scheme was implemented into ABAQUS® version 6.7 
which does not require the input of a fixed stabilization factor mentioned above. The adaptive 
automatic stabilization scheme allows ABAQUS® to automatically increase the damping factor if 
required or reduce the value if the instabilities subside. The result obtained for adaptive automatic 
stabilization and a release tolerance value of 0.2 is plotted in Figure 16. Initially the computed load 
overshot the benchmark result. For increasing propagation, however, the average curve was in good 
agreement with the benchmark result. 
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The results obtained for the coarse meshes shown in Figures 6a (plotted in red) and 6b 
(plotted in blue) are shown in Figure 17. A stabilization factor of 2x1 O’ 1 and a release tolerance of 
0.002 was used for the coarse meshes. One result obtained for the fine mesh with the same 
stabilization factor and a release tolerance of 0.2 (in black) was added to the plot as a reference 
result. Changing the mesh size significantly influenced the magnitude of the saw tooth pattern. 
Larger elements yielded an increased saw tooth in spite of the fact a release tolerance two orders of 
magnitude smaller than the reference result was chosen. 

For the results plotted in Figure 18, contact stabilization was added to the analysis. For all 
combinations of stabilization factors and release tolerances, a saw tooth pattern was observed, 
where the peak values were in good agreement with the benchmark result. The saw tooth curve is 
slightly lower. Decreasing the stabilization factors appeared to cause a slower rate of convergence 
which is either seen by smaller <5/2 for the same number of analysis increments or early tennination 
of the analysis as shown in Table II. Changing the release tolerance also appeared to influence the 
convergence. However, it appeared to have no effect on the overall load/displacement behavior or 
the magnitude of the saw tooth pattern. 

The results obtained for models made of continuum shell elements (shown in Figure 7) are 
shown in Figure 19. A stabilization factor of lxlO" 7 and a release tolerance of 0.2 was used for the 
fine mesh where three continuum shell elements were used over the thickness of one ann as shown 
in Figure 5b (in green). The release tolerance was lowered to a value of 0.002 for the fine mesh 
shown in Figure 7a where only one element was used over the thickness of one ann (in black). The 
initial stiffness of the shell model is slightly reduced. The propagation results, however, are in good 
agreement. It seems that the element type used to model the specimen has no effect on the observed 
saw tooth behavior. A stabilization factor of lxlO' 7 and a release tolerance of 0.2 were used for the 
models with coarser meshes. The results obtained for the models in Figure 7b (plotted in red) and 7c 
(plotted in blue) indicate that changing the mesh size, significantly influenced the magnitude of the 
saw tooth pattern. As observed before, larger elements yielded an increased saw tooth pattern. 

Viscous regularization was added to the analysis to overcome convergence problems. 
Convergence could not be achieved over a wide range of viscosity coefficients when a release 
tolerance value of 0.2 was used as suggested in reference [19]. Subsequently, the release tolerance 
value was increased. The results where convergence was achieved are plotted in Figure 20. For all 
combinations of the viscosity coefficient and release tolerance, a saw tooth pattern was obtained, 
where the peak values were in good agreement with the benchmark result. The average results are 
somewhat lower than the benchmark result. Compared to results obtained from analyses with global 
and contact stabilization, the results obtained with viscous regularization appear to have a better rate 
of convergence since a higher opening displacement ( 5 / 2 = 1.48 mm) was applied during the 
analysis for the same number of total increments (1000). Decreasing the viscosity coefficient 
appeared to cause a slower rate of convergence which was seen by smaller 5/2 values for the same 
number of analysis increments as visible in the plots. Lowering the release tolerance also appeared 
to influence the convergence which was either seen by smaller 5/2 for the same number of analysis 
increments as visible in the plots or early tennination of the analysis as shown in Table II. Changing 
the release tolerance, however, appeared to have no effect on the overall load/displacement behavior 
or the magnitude of the saw tooth pattern. 

In summary, good agreement between analysis results and the benchmark could be achieved 
for different release tolerance values in combination with global or contact stabilization or viscous 
regularization. Selecting the appropriate input parameters, however, was not straightforward and 
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often required several iterations where the parameters had to be changed. All results had a saw 
tooth pattern which appears to depend on the mesh size at the front. 


5.2.2 Computed delamination lengths for DCB specimen with unidirectional layup 

An alternate way to plot the benchmark is shown in Figures 2 1 and 22 where the 
delamination length a is plotted versus the applied opening displacement 8/2 (Figure 21) and the 
computed load P (Figure 22). This way of presenting the results is shown since it may be of 
advantage for large structures where local delamination propagation may have little effect on the 
global stiffness of the structure and may therefore not be visible in a global load/displacement plot. 

For the examples plotted in Figure 23 and 24, global stabilization was used in the analysis. 
For a stabilization factor of 2xl0" 6 (in green), a saw tooth pattern was observed when the 
delamination length a was plotted versus the applied opening displacement 8 / 2 (Figure 23). The 
average, however, was in good agreement with the benchmark result. For a stabilization factor of 
2xl0" 7 (in red), the average appeared to be in better agreement with the benchmark result until 
tennination after 550 increments due to convergence problems. For a stabilization factor of 2x10 s 
(in black for a release tolerance of 0.2), the computed results were on top of the previous result. 
When the delamination length a is plotted versus the computed load P, the saw tooth pattern is more 
pronounced as shown in Figure 24 for the same stabilization factors as above. The average results 
are in good agreement with the benchmark result. 

5.2.3 Computed total strain energy and damping energies for DCB specimen with 
unidirectional layup 

The total strain energy U of the DCB is calculated from the external load P and the crack 
opening displacement 8 as shown in Figure 3 such that 


U = 


P - 5 
2 


( 6 ) 


The calculation is illustrated in Figure 25 using the load/displacement benchmark curve 
discussed above. The areas under the load/displacement curve correspond to the energies U a 
calculated for one ann of the DCB specimen for different loads P and at different delamination 
lengths a. The total strain energies ALLSE obtained from ABAQUS® are plotted in Figure 26 versus 
the applied opening displacement 812 for models of the DCB specimen with different delamination 
length a. The quadratic relationship between the total strain energy and the opening displacement 

U = — and P = C- 6^^-^ (7) 

2 2 


is clearly visible for constant compliance C (constant delamination length a). For comparison, the 
total strain energies U for a model with a = 40 mm were calculated using the applied opening 
displacement and the computed load (equation 7). The results were included in Figure 26 (solid 
brown triangles) and show an excellent agreement with the curve fit through internally computed 
results from ABAQUS® (ALLSE). 

It was assumed that the energies calculated for the critical load/displacement curve can be 
used as a benchmark with respect to the total strain energy. For the delamination propagation, 
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therefore, the results obtained from the model of a DCB specimen with an initially straight 
delamination of a = 30 mm length should follow the benchmark path (in red) in Figure 26. 

In the VCCT for ABAQUS® manual, it is suggested to monitor the energy absorbed by 
damping: ALLSD for contact or global stabilization, ALLVD for viscous damping [19]. The amount 
of damping energy in the models is compared to the total strain energy in the model (ALLSE). 
Ideally, the value of the damping energy should be a small fraction of the total energy. In Figures 27 
to 32, the computed damping energies and the total strain energy in the model of the DCB specimen 
are plotted versus the applied crack tip opening {512) for different input parameters which 
correspond to the results shown in Figures 15 to 20. 

For the results plotted in Figure 27, global stabilization was added to the analysis and the 
analysis results shown correspond to the load/displacement results shown in Figure 15. For a 
stabilization factor of 2xl0" 5 , the calculated total strain (plotted in blue) exceeds the benchmark 
result. For the other stabilization factors of 2xl0" 6 , 2x1 O' 7 , 2x1 0 s (plotted in green, red and black) 
the calculated total strain energies plotted are almost identical and fall slightly below the benchmark 
result. A saw tooth pattern is observed for all the results. For different input parameters, 
significantly different stabilization energies were computed. Below the critical point, the 
stabilization energy was basically zero. Once delamination propagation starts, the stabilization 
energy was required to avoid numerical problems. The lowest stabilization energies were observed 
for stabilization factors of 2xl0" 6 , 2xl0" 7 , 2x1 0 s in combination with a release tolerance of 0.2 
(plotted in green, red and black). The results were almost identical and reached about 20% of the 
total strain energy in the model. Lowering the release tolerance to 0.02 (plotted in light blue) and 
0.002 (plotted in violet) for a stabilization factor of 2x10' appears to increase the stabilization 
energy to about 25% in the example shown. The results obtained for a stabilization factor of 2x10° 
and a release tolerance of 0.2 lie in the middle. 

The result obtained for adaptive automatic stabilization and a release tolerance value of 0.2 
is plotted in Figure 28. The analysis results shown correspond to the load/displacement results 
shown in Figure 16. The calculated total strain energy, for which a saw tooth pattern is observed, 
falls slightly below the benchmark result. For applied displacements below the critical point, the 
stabilization energy was basically zero. Once delamination propagation starts, the stabilization 
energy was required to avoid numerical problems and reached about 25% of the total strain energy 
in the model. 

The results obtained for the coarse meshes shown in Figures 6a (plotted in red) and 6b 
(plotted in blue) are shown in Figure 29. The analysis results shown correspond to the 
load/displacement results shown in Figure 17. A stabilization factor of 2x1 O' 8 and a release 
tolerance of 0.002 was used for the coarse meshes. One result obtained for the fine mesh with the 
same stabilization factor and a release tolerance of 0.2 (in black) was added to the plot as a reference 
result. Changing the mesh size significantly influenced the magnitude of the saw tooth pattern, the 
peak values, however, were in good agreement with the benchmark total strain energy. Larger 
elements yielded an increased saw tooth pattern in spite of the fact a release tolerance two orders of 
magnitude smaller than the reference result was chosen. The lowest stabilization energy (about 20% 
of the total strain energy) was observed for the reference result with a stabilization factor of 2xl0" 8 
in combination with a release tolerance of 0.2. Lowering the release tolerance to 0.002 and 
increasing the element size appears to increase the stabilization energy to more than 50% in the 
example shown. 

For the results plotted in Figure 30, contact stabilization was added to the analysis. These 
results correspond to the load/displacement results shown in Figure 18. For all combinations of 
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stabilization factors and release tolerances, a saw tooth pattern was observed, where the peak values 
were in good agreement with the benchmark total strain energy. The lowest stabilization energies 
were observed for stabilization factors of lxlO" 5 , lxlO" 6 , lxlO’ 7 in combination with a release 
tolerance of 0.2 (plotted in red, blue and light blue). The stabilization energy results were almost 
identical and reached more than 20% of the total strain energy in the model. Lowering the release 
tolerance to 0.02 (plotted in orange) and 0.002 (plotted in green) for a stabilization factor of lxlO" 7 
appears to increase the stabilization energy to more than 30% in the example shown. The 
stabilization energy followed the same path for a release tolerance of 0.002 for stabilization factors 
of lxlO’ 7 and lxlO’ 3 . 

The results obtained for models made of continuum shell elements SC8R (shown in 
Figure 7) are shown in Figure 31. The analysis results shown correspond to the load/displacement 
results shown in Figure 19. A stabilization factor of lxlO" 7 and a release tolerance of 0.2 were used 
for the fine mesh where three continuum shell elements were used over the thickness of one ann as 
shown in Figure 5b (in green). The release tolerance was lowered to a value of 0.002 for the fine 
mesh shown in Figure 7a where only one element was used over the thickness of one arm. The 
calculated total strain energies plotted are almost identical and fall slightly below the benchmark 
result. Changing the mesh size for a stabilization factor of lxlO’ 7 and a release tolerance of 0.2 
significantly influenced the magnitude of the saw tooth pattern, the peak values, however, were in 
good agreement with the benchmark total strain energy. Larger elements yielded an increased saw 
tooth pattern. The lowest stabilization energy (about 20% of the total strain energy) was observed 
for the fine mesh shown in Figure 5b, a stabilization factor of lxlO" 7 and a release tolerance value of 
0.2 (plotted in green). Larger elements yielded an increased saw tooth pattern (plotted in red and 
blue) but also increased the stabilization energy required to more than 35%. Lowering the release 
tolerance to a values of 0.002 (plotted in black) also lead to an increase in stabilization energy. 

For all combinations of the viscosity coefficient and release tolerance, a saw tooth pattern 
was obtained, where the peak values were in good agreement with the benchmark total strain energy 
as shown in Figure 32. The analysis results shown correspond to the load/displacement results 
shown in Figure 20. The lowest stabilization energies were observed for stabilization factors of 
lxlO" 4 , lxlO' 5 , in combination with a release tolerance of 0.5 (plotted in red and green). The results 
were almost identical and reached only about 5% of the total strain energy in the model. Lowering 
the release tolerance to 0.3 for viscosity coefficients of lxlO" 4 , lxlO" 3 (plotted in blue and black) 
appears to increase the stabilization energy to more than 15% in the example shown. 

In summary, good agreement between analysis results and the total strain energy benchmark 
could be achieved for different release tolerance values in combination with global or contact 
stabilization or viscous regularization. All results had a saw tooth pattern and the magnitude of 
which appears to depend on the mesh size at the delamination front. Larger elements yielded an 
increased saw tooth pattern. Stabilization energies of about 20%-25% of the total stain energy were 
observed when release tolerance values of 0.5 and 0.2 were used. Lowering the release tolerance to 
values of 0.02 or 0.002 resulted in an increase in stabilization energy. The lowest stabilization 
energies (about 5% of the total strain energy in the model) were observed for viscosity coefficients 
of lxlO" 4 , lxlO" 3 , in combination with a release tolerance of 0.5. In spite of the variations in 
stabilization energy, all the load/displacement results were in good agreement with the benchmark 
as shown earlier in Figures 15 to 20. It is therefore uncertain if the amount of stabilization energy 
absorbed can be used as a measure to detennine the quality of the analysis results. 
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5.2.4 Computed delamination front shape 

Besides matching the load displacement behavior of benchmark results, a delamination 
propagation analysis should also yield a delamination front shape that is representative of the 
actual failure. An example of delamination front shapes observed by opening a tested DCB 
specimen are shown in Figure 33a [43]. From the initial straight delamination front which is fonned 
by the edge of the Teflon insert, the delamination develops into a curved thumbnail shaped front. 
The front remains thumbnail shaped if the test is continued and the delamination continues to grow. 
Delamination propagation computed using the model with a uniform mesh across the width 
(Figure 5b) is shown in Figure 33b at the end of the analysis after 1000 increments. Plotted on the 
bottom surface (defined in Figure 5b) are the contours of the bond state, where the delaminated 
section appears in red and the intact (bonded) section in blue. The transition between the colors 
indicates the location of the delamination front. The initial straight front was included for 
clarification. The first propagation was observed in the center of the specimen as expected from the 
distribution of the energy release rate (Figure 11) and the failure index (Figure 12). The front 
propagated across the width of the specimen until a new straight front was reached. Subsequently, 
the propagation starts again in the center. During the analysis, the front never developed into the 
expected curved thumbnail front, and the analysis tenninated with a straight front as shown in 
Figure 33b. This result is somewhat unsatisfactory but may be explained by the fact that the failure 
index in this particular example is nearly constant across about 80% of the width of the specimen as 
shown in Figure 12. An even finer mesh may be required to capture the lagging propagation near 
the edge. 


5.3 Creating a Benchmark Solution for a DCB specimen with multi-directional layup 

The analysis outlined in Section 4.1 were repeated for a DCB specimen with multi- 
directional layup. First, the mode I strain energy release rate values were computed which are 
plotted versus the normalized width, y/B, of the specimen as shown in Figure 34. The results were 
obtained from models shown in Figure 8 for ten different delamination lengths a. An opening 
displacement d/2=1.0 mm was applied to each arm of the model. Qualitatively, the mode I strain 
energy release rate is fairly constant in the center part of the specimen and drops progressively 
towards the edges. Compared to the DCB with unidirectional layup, the constant center section is 
smaller, and the drop towards the edges occurs earlier for specimens with the multi-directional 
layup. These effects are caused by increased anticlastic bending in the more compliant specimens 
with the multi-directional layup [24, 25], Computed mode I strain energy release rates decreased 
with increasing delamination length a. 

The failure index Gt/G c was computed next, based on a mode I fracture toughness 
Gj c = 340.5 J/m 2 for C12K/R6376 (see Figure 2b). The failure index was plotted versus the 
nonnalized width, y/B, of the specimen as shown in Figure 35. For all delamination lengths 
modeled, except for <7=40 mm, the failure index in the center of the specimen (y/B= 0) is above unity 
(G t /G c >\). 

For all delamination lengths modeled, the reaction loads P at the location of the applied 
displacement were calculated and plotted versus the applied opening displacement d/2 as shown in 
Figure 36. The critical load P cr it and critical displacement <5 ni ,/2 were calculated for each 
delamination length modeled using equation (5), and the results were included in the 
load/displacement plots as shown in Figure 37 (solid red circles). As before, it was assumed that the 
critical load/displacement results can be used as a benchmark. For the delamination propagation, 
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therefore, the load/displacement results obtained from the model of a DCB specimen with an 
initially straight delamination of a = 30 mm length should correspond to the critical 
load/displacement path (solid red line) in Figure 37. 


5.4 Delamination Propagation in a DCB specimen with multi-directional layup 
5.4.1 Computed load/displacement behavior for different input parameters 

The propagation analysis was performed in two steps using the model shown in Figure 8 for 
a delamination length 3 1 mm. In the first step, a prescribed displacement {812= 0.7 mm) was applied 
in two increments. Dividing the first step into just two increments was possible, since the load- 
displacement behavior of the specimen up to failure was linear as shown in Figure 37. In the second 
step, the total prescribed displacement was increased {812= 2.8 mm). Automatic incrementation was 
used with a small increment size at the beginning (10 4 of the total increment) and a very small 
minimum allowed increment (10’ of the total increment) to reduce the risk of numerical instability 
and early tennination of the analysis. The analysis was limited to 1000 increments. 

In Figure 38, the computed resultant force (load P) at the tip of the DCB specimen is plotted 
versus the applied crack tip opening {812) for different input parameters. For the results shown, the 
analysis tenninated when the 1000 increment limit set for the analysis was reached. To overcome 
the convergence problems, global stabilization, contact stabilization and viscous regularization were 
used with input parameters for which good results had been obtained for the analysis of the 
unidirectional DCB specimen. Viscous regularization was added to the analysis with a viscosity 
coefficient of lxlO 0 and a release tolerance of 0.3 (plotted in green). The analysis using contact 
stabilization was performed with a stabilization a factor of lxlO" 7 and a release tolerance of 0.002 
(plotted in blue). For global stabilization, a factor of 2xl0’ 7 was used in combination with a release 
tolerance value of 0.2 (plotted in red). For all different input parameters, which are also listed in 
Table II, converged solutions were obtained, and the analyses reached the predetennined 1000 
increment limit. All results had a saw tooth pattern, where the average values were in good 
agreement with the benchmark result. Compared to results obtained from analyses with global 
stabilization (in red) and contact stabilization (in blue), the results obtained with viscous 
regularization (in green) appear to have a better rate of convergence. For the same number of total 
increments (1000), the analysis continued to a higher opening displacement {8/2= 1.24 mm). 


5.4.2 Computed total strain energy and damping energies for DCB specimen with multi- 
directional layup 

As discussed earlier for the DCB specimen with unidirectional layup, it is assumed that the 
plot of the total strain energy versus the applied crack tip opening displacement {8/ 2) can be used as 
a benchmark result. The benchmark result obtained from ABAQUS® analyses from models of the 
multi-directional DCB specimen is plotted (in grey) in Figure 39. For the delamination propagation, 
the results obtained from the model of a DCB specimen with an initially straight delamination of 
< 7=3 1 mm length are expected to follow the path of the benchmark. 

In Figure 39, the total strain energy in the model of the DCB specimen and the computed 
damping energies are plotted versus the applied crack tip opening {8/ 2) for different input 
parameters which correspond to the results shown in Figure 38. For all the results plotted, a saw 
tooth pattern is observed. The calculated total strain energies are in good agreement with the 
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benchmark result. For a viscosity coefficient of lxlO" 5 and a release tolerance of 0.3 (plotted in 
green), the total strain energy slightly exceeds the benchmark result. For the different input 
parameters, significantly different stabilization energies were computed. Below the critical point, 
the stabilization energy was basically zero for all results. Once delamination propagation starts, the 
stabilization energy was required to avoid numerical problems. The lowest stabilization energies 
were observed for a viscosity coefficient of lxl O’ 5 and a release tolerance of 0.3 (plotted in green) 
and reached about 1 1% of the total strain energy in the model. Lowering the release tolerance to 0.2 
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(plotted in red) for a global stabilization factor of 2x10" and 0.002 (plotted in blue) for a contact 
stabilization factor of lxlO" 7 appears to increase the stabilization energy to about 20% in the 
example shown. In spite of the variations in stabilization energy, all the load/displacement results 
were in good agreement with the benchmark. It is therefore uncertain if the amount of stabilization 
energy absorbed can be used as a measure to detennine the quality of the analysis results. 


5.4.3 Computed delamination front shape for a DCB specimen with multi-directional layup 

An initial straight delamination front develops into a curved thumbnail shaped front which 
was shown for an opened tested DCB specimen in Figure 33a. The front remains thumbnail shaped 
if the test is continued and the delamination continues to grow. The thumbnail shaped front is 
caused by the anticlastic bending of the anns which is more prevalent in the more compliant arms of 
a multi-directional specimen [44], Delamination propagation computed using the model with a 
uniform mesh across the width (Figure 8) is shown in Figure 40 at the end of the analysis after 1000 
increments. Plotted on the bottom surface (defined in Figure 8) are the contours of the bond state, 
where the delaminated section appears in red, and the intact (bonded) section in blue. The transition 
between the colors indicates the location of the delamination front. The initial straight delamination 
front was included for clarification. The first propagation was observed in the center of the 
specimen as expected from the distribution of the energy release rate (Figure 34) and the failure 
index (Figure 35). During the analysis, the front developed into the expected curved thumbnail 
front. Compared to results obtained from analyses with contact stabilization (Figure 40a), the results 
obtained with viscous regularization (Figure 40b) appear to have a better rate of convergence since 
for the same limit of 1000 increments the front grew further into the specimen. Compared to the 
straight fronts obtained from the models of the unidirectional DCB specimen (shown in Figure 33b), 
the current results are encouraging. It remains however somewhat unclear what degree of mesh 
refinement is required to accurately capture the delamination front shape. 


5.5 Creating a Benchmark Solution for a SLB specimens 

The computed total strain energy release rate values were plotted versus the nonnalized 
width, y/B, of the SLB specimen as shown in Figure 41. The results were obtained from 
geometrically nonlinear analysis of models shown in Figure 9 for twelve different delamination 
lengths a. An arbitrary center deflection w= 2.8 mm was applied as shown in Figure 3. Qualitatively, 
the total energy release rate is fairly constant in the center part of the specimen and drops towards 
the edges. Peaks in the distribution are observed at the edges. Computed total strain energy release 
rates decreased with increasing delamination length a. 

The sum of the shear components Gs = Gn+Gni and the mixed-mode ratio Gs /GY were also 
calculated for each nodal point along the delamination front across the width of the specimen. The 
mixed-mode ratio Gs/Gt was plotted versus the normalized width, y/B, of the specimen as shown in 
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Figure 42. Qualitatively, the mixed-mode ratio is fairly constant in the center part of the specimen 
progressively increasing towards the edges. Using the mixed-mode failure criterion for 
C12K/R6376 (see Figure 2b), the failure index Gt/G c was computed for each node along the 
delamination front and plotted versus the nonnalized width, y/B, of the specimen as shown in 
Figure 43. For the center deflection applied, the failure index Gj/G, in the center is well below 
unity. The failure index is almost constant in the center of the specimen, drops towards the edges 
and increases again in the immediate vicinity of the edge. To reach Gj/G c = 1 in the center of the 
specimen (y/B= 0), a critical center deflection, w cnh and corresponding critical load P cr it, were 
calculated using equation (5) for all delamination lengths modeled. 

For all delamination lengths modeled, the reaction load P at the location of the applied 
deflection were calculated and plotted versus the applied center deflection, w, as shown in 
Figure 44. The calculated critical center deflection, w cr u, and corresponding critical load values, P cr it, 
were added to the plots as shown in Figure 45 (solid red circles). The results indicated that, with 
increasing delamination length, less load is required to extend the delamination. At the same time 
also, the values of the critical center deflection decreased. This means that the SLB specimen 
exhibits unstable delamination propagation under load as well as displacement control (dashed red 
line). From these critical load/displacement results, a benchmark solution can be created. To define 
the benchmark, it is assumed that prescribed center deflections are applied in the analysis instead of 
nodal point loads P to minimize problems with numerical stability of the analysis caused by the 
unstable propagation. Once the critical center deflection is reached and delamination propagation 
starts, the applied displacement must be held constant over several increments while the 
delamination front is advanced during these increments. Once the stable path is reached, the applied 
center deflection is increased again incrementally. For the simulated delamination propagation, 
therefore, the load/displacement results obtained from the model of a SLB specimen with an 
initially straight delamination length of a=34 mm should correspond to the benchmark 
load/displacement path (solid red line) as shown in Figure 45. 


5.6 Delamination Propagation in a SLB Specimen using VCCT for ABAQUS® 

5.6.1 Computed load/displacement behavior for different input parameters 

The propagation analysis was perfonned in two steps using the models shown in Figures 8 
and 9. In the first step, a central deflection ( w= 3.1 mm) was applied in two increments which 
equaled nearly the critical center deflection (w cr u = 3.23 mm) detennined in the analysis above. 
Dividing the first step into just two increments was possible, since the load-displacement behavior 
of the specimen up to failure was linear as shown in Figure 45. In the second step, the total 
prescribed displacement was increased ( w= 5.0 mm). Automatic time incrementation was used with 
a small initial time increment size (10 3 ) and a very small minimum allowed time increment (10 17 ) 
to reduce the risk of numerical instability and early termination of the analysis. The analysis was 
limited to 1000 increments. 

In Figures 46 to 50, the computed resultant force (load P ) at the center of the SLB specimen 
is plotted versus the center deflection (w) for different input parameters which are listed in Table II. 
The analysis terminated before the total prescribed center deflection was applied. For the results 
shown, the analysis terminated when the 1000 increment limit set for the analysis was reached. 
Several analyses tenninated early because of convergence problems. The results computed when 
global stabilization was used are plotted in Figure 46. For a stabilization factor of 2x1 0’ 5 , the load 
increased suddenly at the beginning of the second load step (plotted in blue). Then, the load 
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continued to increase on a path with the same stiffness as the benchmark but offset to higher loads. 
The load continued to increase until a point was reached where delamination propagation started 
and the load decreased. The analysis was stopped by the user. For a stabilization factor of 2xl0" 6 (in 
green), the delamination propagation started at the critical center deflection. In the beginning, the 
load/displacement path followed the constant deflection branch of the benchmark result very well. 
At the transition between the constant deflection branch and the stable propagation branch of the 
benchmark result, the applied center deflection was about 2% higher compared to the benchmark. 
For the stable path, a saw tooth pattern was observed but the minimum is in good agreement with 
the benchmark result. 

An adaptive automatic stabilization scheme was recently added to ABAQUS®/Standard in 
version 6.7. The adaptive automatic stabilization scheme does not require the input of a fixed 
stabilization factor mentioned above. The results obtained are plotted in Figure 47. For the default 
setting, the load increased suddenly at the beginning of the second load step (plotted in blue). Then, 
the load continued to increase on a path with the same stiffness as the benchmark but offset to 
higher loads. The load continued to increase until a point was reached where delamination 
propagation started and the load decreased. The analysis was stopped by the user. To obtain 
converged results, the stabilization factor at the beginning of the analysis was detennined by the 
user and automatic stabilization adjusted the settings in the following increments. For an initial 
stabilization factor of 2xl0" 6 (in red), the delamination propagation started at the critical center 
deflection. In the beginning, the load/displacement path followed the constant deflection branch of 
the benchmark result very well but overshot and did not follow the stable path. When the initial 
stabilization factor was changed to 2xl0" 7 (in light blue) and 2xl0" 8 (in black), the same path was 
followed as before but the analysis tenninated early. Increasing the release tolerance to 0.5 and 0.9 
for selected initial stabilization factors 2x1 O’ 7 (in orange) and 2xl0" 6 (in green) did not lead to a 
converged solution. Further improvement of the automatic stabilization scheme is required before it 
can be used reliably. 

The results computed when contact stabilization was used are plotted in Figure 48. For a 
small stabilization factor (lxlO" 6 ) and a release tolerance (0.2) suggested in the handbook [19], the 
load dropped and delamination propagation started prior to reaching the critical point of the 
benchmark solution (plotted in blue). The load/displacement path then ran parallel to the constant 
deflection branch of the benchmark result but the analysis tenninated early due to convergence 
problems. The stabilization factor and release tolerance had to be increased to avoid premature 
termination of the analysis. For a stabilization factor of lxlCT and release tolerance of 0.5 (in 
green), the load dropped at the critical point of the benchmark solution. First, the center deflection 
kept increasing with decreasing load. Later, the load/displacement path ran parallel to the constant 
deflection branch of the benchmark result. At the transition between the constant deflection branch 
and the stable propagation branch of the benchmark result, the applied center deflection was about 
2% higher compared to the benchmark. For the stable path, a saw tooth pattern was observed where 
the average results were in good agreement with the benchmark result. The difference between the 
maximum and minimum values was much smaller than in the case where global stabilization was 
used (see Figure 46). The best results compared to the benchmark were obtained for even higher 
values of the stabilization factors of lxlO" 4 and a release tolerance of 0.5 (in red). 

When viscous regularization was used to help overcome convergence issues, a value of 0.2 
was used initially for the release tolerance as suggested in the handbook [19]. Convergence could 
not be achieved which led to an increase in the release tolerance. The results are plotted in 
Figure 49. For a small viscosity coefficient of lxlO" 4 and a release tolerance of 0.5 (in blue), the 
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load dropped at the critical point, but the center deflection kept increasing with decreasing load. 
Then, the analysis tenninated early due to convergence problems. For an increased viscosity 
coefficient of 1 x 1 0" 2 and a release tolerance of 0.5 (in red), the load dropped at the critical point and 
the load/displacement path started following the constant deflection branch of the benchmark result, 
but the analysis tenninated early due to convergence problems. The viscosity coefficient and release 
tolerance had to be increased further to avoid premature tennination of the analysis. For a viscosity 
coefficient of lxlO" 1 and a release tolerance of 0.9 (in green), the load dropped at the critical point. 
First, the center deflection kept increasing with decreasing load. Later, the load/displacement path 
ran parallel to the constant deflection branch of the benchmark result. At the transition between the 
constant deflection branch and the stable propagation branch of the benchmark result, the applied 
center deflection is about 2.5% higher compared to the benchmark. For the stable path, a saw tooth 
pattern is observed where the average results are in good agreement with the benchmark result. The 
difference between the maximum and minimum values is much smaller compared to the cases 
where global or contact stabilization was used. 

In Figure 50, the computed resultant force (load P ) in the center of the SLB specimen is 
plotted versus the applied center deflection (w) for the model of the SLB specimen shown in 
Figure 10. For the results shown, the analysis terminated when the 1000 increment limit set for the 
analysis was reached. To overcome the convergence problems, global stabilization, contact 
stabilization and viscous regularization were used with input parameters for which good results had 
been obtained previously (FE model shown in Figure 9). For a viscosity coefficient of lxlO" 1 and a 
release tolerance of 0.9 (in green), the load dropped at the critical point. First, the center deflection 
kept increasing with decreasing load. Later, the computed load/displacement path ran parallel to the 
constant deflection branch of the benchmark result. At the transition between the constant deflection 
branch and the stable propagation branch of the benchmark result, the applied center deflection was 
about 2.5% higher compared to the benchmark. For contact stabilization, a factor of lxlO" 4 and a 
release tolerance of 0.5 were chosen (plotted in blue). The computed load started decreasing at the 
critical point of the benchmark solution. First, the center deflection kept increasing with decreasing 
load. Later, the computed path ran parallel to the constant deflection branch of the benchmark result. 
At the transition between the constant deflection branch and the stable propagation branch of the 
benchmark result, the applied center deflection was about 2% higher compared to the benchmark. 
For a stabilization factor of 2xl0" 6 (in red), the delamination propagation started at the critical center 
deflection. In the beginning, the load/displacement path followed the constant deflection branch of 
the benchmark result very well. At the transition between the constant deflection branch and the 
stable propagation branch of the benchmark result, the applied center deflection was about 2% 
higher compared to the benchmark. For all results, a saw tooth pattern was observed along the stable 
path, the average results however, were in good agreement with the benchmark result. 

In summary, good agreement between analysis results and the benchmark could be achieved 
for different release tolerance values in combination with global or contact stabilization or viscous 
regularization. Selecting the appropriate input parameters, however, was not straightforward and 
often required several iterations in which the parameters had to be changed. 


5.6.2 Computed delamination lengths for SLB specimen with multi-directional layup 

An alternate way to plot the benchmark is shown in Figures 5 1 and 52 where the 
delamination length a is plotted versus the applied center deflection w (Figure 51) and the computed 
load P (Figure 52). This way of presenting the results is shown, since it may be of advantage for 


19 



large structures where local delamination propagation may have little effect on the global stiffness 
of the structure and may therefore not be visible in a global load/displacement plot. 

For the examples plotted in Figures 53 and 54, global stabilization and contact stabilization 
was used in the analysis. For a global stabilization factor of 2xl0" 6 (in green), the delamination 
propagation started at the critical center deflection as shown in Figure 53. The plotted results 
followed the benchmark result very well including the transition between the constant deflection 
branch and the stable propagation branch of the benchmark result. For the stable path, a saw tooth 
pattern was observed where the maximum was in good agreement with the benchmark result. For a 
contact stabilization factor of lxlO" 4 and a release tolerance of 0.5 (in red) the results also followed 
the benchmark result very well. When the delamination length a is plotted versus the computed load 
P, the saw tooth pattern becomes more visible as shown in Figure 54 for the same stabilization 
factors as above. The increased saw tooth pattern corresponds to the pattern observed in the 
load/displacement plot shown in Figures 46 and 48. In spite of the saw tooth pattern observed, the 
average results are in good agreement with the benchmark result. 


5.6.3 Computed total strain energy and damping energies for SLB specimen with multi- 
directional layup 

As discussed earlier for the DCB specimen, it is assumed that the total strain energy versus 
applied center deflection plot can also be used as a benchmark result. For models of the SLB 
specimen with different delamination length a, the total strain energies ALLSE obtained from 
ABAQUS® are plotted versus the applied center deflection, w, as shown in Figure 55. The 
calculated total strain energy values corresponding to the critical center deflection, w C nt, and critical 
load values, P a -it, at delamination onset (shown in Figure 45) were added to the plots (solid red 
circles). The results indicated that, with increasing delamination length, the total strain energy in the 
system decreases. Simultaneously, the values of the critical center deflection also decrease. This 
means that the SLB specimen exhibits unstable delamination propagation under load as well as 
displacement control. From these critical energy/displacement results, a benchmark solution can be 
created. To define the benchmark, it is assumed that prescribed center deflections are applied in the 
analysis instead of nodal point loads P to minimize problems with numerical stability of the analysis 
caused by the unstable propagation. Once the critical center deflection is reached and delamination 
propagation starts, the applied displacement must be held constant over several increments while the 
delamination front is advanced during these increments. Once the stable path is reached, the applied 
center deflection is increased again incrementally. For the simulated delamination propagation, 
therefore, the total strain energy/displacement results obtained from the model of a SLB specimen 
with an initially straight delamination length of a = 34 mm should correspond to the critical total 
strain energy/displacement path (solid red line) as shown in Figure 55. 

In Figure 56, the total strain energy in the model of the SLB specimen and the computed 
global stabilization energies are plotted versus the applied center deflection (w) for different input 
parameters which correspond to the results shown in Figure 46. For the two results plotted, a saw 
tooth pattern is observed. For a stabilization factor of 2x10° (in blue), the total strain energy first 
followed the benchmark result but continued to increase once the critical point was reached. The 
energy increased until a point was reached where delamination propagation started and the energy 
started decreasing. The analysis was stopped by the user. For a stabilization factor of 2xl0" 6 (in 
green), the delamination propagation started at the critical center deflection. In the beginning, the 
total strain energy/displacement path followed the constant deflection branch of the benchmark 
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result very well. At the transition between the constant deflection branch and the stable propagation 
branch of the benchmark result, the applied center deflection was about 2% higher compared to the 
benchmark. For the stable path, a saw tooth pattern was observed but the average energy is in good 
agreement with the benchmark result. The stabilization energy was basically zero for both results 
before delamination onset. Once delamination propagation starts, the stabilization energy was 
required to avoid numerical problems. For a stabilization factor of 2x1 O’ 6 (in green), the stabilization 
energy reached about 40% of the total strain energy in the model, the load/displacement results 
plotted in Figure 46 however were in good agreement with the benchmark. 

In Figure 57, the total strain energy in the model of the SLB specimen and the computed 
stabilization energies are plotted versus the applied center deflection (w) for different input 
parameters which correspond to the results for automatic stabilization shown in Figure 47. For the 
default factor (in blue), the total strain energy first followed the benchmark result but continued to 
increase once the critical point was reached. The energy increased until a point was reached where 
delamination propagation started and the energy started decreasing. The analysis was stopped by the 
user. None of the other settings yielded a converged solution. The results up to the tennination of 
the analysis are included for completeness. The stabilization energy was basically zero before 
delamination onset. Once delamination propagation starts, the stabilization energy was required to 
avoid numerical problems. Since none of the analyses reached the stable branch of the benchmark 
before termination, it is unclear how much stabilization energy would have been required. As 
mentioned above, further improvement of the automatic stabilization scheme is required before it 
can be used reliably. 

The results computed when contact stabilization was used are plotted in Figure 58. The 
results correspond to the load/displacement plots shown in Figure 48. For a small stabilization factor 
(lxlO" 6 ) and a release tolerance (0.2) suggested in the handbook [19], delamination propagation 
started prior to reaching the critical point of the benchmark solution (plotted in blue), and the energy 
dropped. The path then ran parallel to the constant deflection branch of the benchmark result but the 
analysis terminated early due to convergence problems. The stabilization factor and release 
tolerance were increased to avoid premature tennination of the analysis. For a stabilization factor of 
lxlO" 3 and release tolerance of 0.5 (in green), the energy started decreasing at the critical point of 
the benchmark solution. First, the center deflection kept increasing with decreasing total strain 
energy. Later, the computed path ran parallel to the constant deflection branch of the benchmark 
result. At the transition between the constant deflection branch and the stable propagation branch of 
the benchmark result, the applied center deflection was about 2% higher compared to the 
benchmark. For the stable path, a saw tooth pattern was observed where the average results were in 
good agreement with the benchmark result. The best results compared to the benchmark were 
obtained for even higher values of the stabilization factors of lxlO" 4 and a release tolerance of 0.5 
(in red). For both results, the stabilization energy reached about 45% of the total strain energy in the 
model, the load/displacement results plotted in Figure 48 however were in good agreement with the 
benchmark. 

The results computed when viscous regularization was used are plotted in Figure 59. The 
results correspond to the load/displacement plots shown in Figure 49. For a small viscosity 
coefficient of lxlO' 4 and a release tolerance of 0.5 (in blue), the calculated total strain energy started 
decreasing at the critical point, but simultaneously the center deflection kept increasing with 
decreasing energy. Then, the analysis was tenninated early due to convergence problems. For an 
increased viscosity coefficient of lxlO" 2 and a release tolerance of 0.5 (in red), the total strain energy 
dropped at the critical point and the computed path started following the constant deflection branch 
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of the benchmark result, but the analysis tenninated early due to convergence problems. The 
viscosity coefficient and release tolerance had to be increased further to avoid premature 
tennination of the analysis. For a viscosity coefficient of lxlO" 1 and a release tolerance of 0.9 (in 
green), the total strain energy dropped at the critical point. First, the center deflection kept 
increasing with decreasing load. Later, the computed path ran parallel to the constant deflection 
branch of the benchmark result. At the transition between the constant deflection branch and the 
stable propagation branch of the benchmark result, the applied center deflection is about 2.5% 
higher compared to the benchmark. For the stable path, a saw tooth pattern was observed where the 
average results are in good agreement with the benchmark result. The stabilization energy reached 
about 34% of the total strain energy in the model, the load/displacement results plotted in Figure 49 
however were in good agreement with the benchmark. 

In Figure 60, the results obtained from the model of the SLB specimen shown in Figure 10 
are plotted versus the applied center deflection, w, for different input parameters. The results 
correspond to the load/displacement plots shown in Figure 50. For a viscosity coefficient of lxlO" 1 
and a release tolerance of 0.9 (in green), the total strain energy dropped at the critical point. First, 
the center deflection kept increasing with decreasing load. Later, the computed path ran parallel to 
the constant deflection branch of the benchmark result. At the transition between the constant 
deflection branch and the stable propagation branch of the benchmark result, the applied center 
deflection is about 2.5% higher compared to the benchmark. For contact stabilization, a factor of 
lxlO" 4 and release tolerance of 0.5 were chosen. The computed total strain energy (in blue) started 
decreasing at the critical point of the benchmark solution. First, the center deflection kept increasing 
with decreasing total strain energy. Later, the computed path ran parallel to the constant deflection 
branch of the benchmark result. At the transition between the constant deflection branch and the 
stable propagation branch of the benchmark result, the applied center deflection was about 2% 
higher compared to the benchmark. For a stabilization factor of 2x1 O’ 6 (in red), the delamination 
propagation started at the critical center deflection. In the beginning, the total strain 
energy/displacement path followed the constant deflection branch of the benchmark result very 
well. At the transition between the constant deflection branch and the stable propagation branch of 
the benchmark result, the applied center deflection was about 2% higher compared to the 
benchmark. For all results, a saw tooth pattern was observed along the stable path, the average 
energy however, was in good agreement with the benchmark result. The stabilization energy before 
delamination onset was basically zero for all results. Once delamination propagation started, the 
stabilization energy was required to avoid numerical problems. The stabilization energy reached 
about 40% of the total strain energy in the model, the load/displacement results plotted in Figure 50 
however were in good agreement with the benchmark. 

In summary, good agreement between analysis results and the total strain energy benchmark 
could be achieved for different release tolerance values in combination with global or contact 
stabilization or viscous regularization. All results had a saw tooth pattern the magnitude of which 
appears to depend on the mesh size at the front. Stabilization energies of about 34%-45% of the 
total stain energy were observed when release tolerance values of 0.9, 0.5 and 0.2 were used. In 
spite of the variations in stabilization energy, all the load/displacement results were in good 
agreement with the benchmark as discussed earlier and shown Figures 46 and 48-50. It is therefore 
uncertain if the amount of stabilization energy absorbed can be used as a measure to detennine the 
quality of the analysis results. 
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5.6.4 Computed delamination front shape for a SLB specimen with multi-directional layup 

Delamination propagation computed using the model with a unifonn mesh across the width 
(Figure 9) is shown in Figures 61 and 62 after 76 and 1000 increments, respectively. Plotted on the 
bottom surface (defined in Figure 9) are the contours of the bond state variable. The bond state 
varies between 0.0 (fully bonded shown in dark blue) and 1.0 (fully disbonded shown in red) [19]. 
The transition between the colors indicates the location of the delamination front. The initial straight 
front was included for clarification. The first propagation is observed near the center and 
corresponds to the maximum in the distribution of the failure index (Figure 43). The front then 
propagated across the width. Further propagation created a curved front where the edges lag behind 
as shown in Figure 61. This result is in good agreement with expectations based on the distribution 
of the failure index shown in Figure 43. Once the delamination propagation extends into the coarser 
part of the mesh, the curved front cannot be represented as well as shown in Figure 62. 
Delamination propagation computed using the model with an extended finer mesh (Figure 10) is 
shown in Figure 63a and b after 476 and 1000 increments, respectively. During propagation through 
the fine section the curved front is clearly visible as shown in Figure 63a. As mentioned above the 
curved front cannot be represented as well once delamination propagation extends into the coarser 
part of the mesh as shown in Figure 63b. C-scans or x-ray photographs of tested specimens were 
not available for comparison. 


6. SUMMARY AND CONCLUSIONS 

An approach for assessing the delamination propagation simulation capabilities in 
commercial finite element codes is presented and demonstrated for the commercial finite element 
code ABAQUS' with focus on their implementation of the Virtual Crack Closure Technique 
(VCCT). For this investigation, the Double Cantilever Beam (DCB) specimen with a unidirectional 
and a multi-directional layup and the Single Leg Bending (SLB) specimen with a multi-directional 
layup were chosen for full three-dimensional finite element simulations. First, critical 
load/displacement results were defined for delamination onset which were used subsequently as 
benchmarks. Second, starting from an initially straight front, the delamination was allowed to 
propagate based on the algorithms implemented into VCCT for ABAQUS 8 . VCCT control 
parameters were varied to study the effect on the computed load-displacement behavior during 
propagation. It was assumed that for good results, the computed load-displacement relationship 
should correspond to the benchmark results. Third, as a qualitative assessment, the shape of the 
computed delamination fronts were also compared to photographs of failed specimens. 

The results showed the following: 

• Good agreement between the load-displacement relationship obtained from the 
propagation analysis results and the benchmark results could be achieved by selecting 
the appropriate input parameters. However, selecting the appropriate VCCT input 
parameters such as release tolerance, global or contact stabilization and viscous 
regularization, was not straightforward and often required an iterative procedure. 

• The default setting for global stabilization yielded unsatisfactory results although the 
analysis converged. 

• Analysis with automatic global stabilization yielded unsatisfactory results. 

• Best results were obtained when contact stabilization and viscous regularization were 
used. 

• Accurately computing the delamination front shape requires fine meshes. 
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• A coarse mesh could lead to an increased saw tooth behavior during propagation 
analysis. 

• Good agreement between analysis results and the total strain energy benchmark could be 
achieved for different release tolerance values in combination with global or contact 
stabilization or viscous regularization. 

• Load/displacement results were in good agreement with the benchmark for cases where 
the stabilization energies varied between 5% and 25% for the analysis of the DCB 
specimen and between 34% and 45% for the analysis of the SLB specimen. It is 
therefore uncertain if the amount of stabilization energy absorbed can be used as a 
measure to determine the quality of the analysis results. 

Overall, the results are promising. In a real case scenario, however, where the results are 
unknown, obtaining the right solution will remain challenging. Further studies are required which 
should include different levels of mesh refinement, new stabilization options and the use of 
continuum shell elements to more efficiently model the specimens. Additionally, assessment of the 
propagation capabilities in more complex specimens and on a structural level is required. 
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TABLE 1 . MATERIAL PROPERTIES. 


T 300/1076 Unidirectional Graphite/Epoxy Prepreg 

E\ i = 139.4 GPa 

£22 = 10.16 GPa 

£33 = 10.16 GPa 

V 12 = 0.30 

vi 3 = 0.30 

V 23 = 0.436 

Gi2 = 4 . 6 GPa 

Gl3 = 4.6 GPa 

G23 = 3.54 GPa 

C 12 K/R 6376 Unidirectional Graphite/Epoxy Prepreg 

Ell = 146.9 GPa 

£22 = 10.6 GPa 

£33 = 10.6 GPa 

V12 = 0.33 

vi3 = 0.33 

V 23 = 0.33 

G12 = 5.45 GPa 

G13 = 5.45 GPa 

G23 = 3.99 GPa 


The material properties are given with reference to the ply coordinate axes where index 1 1 denotes the ply 
principal axis that coincides with the direction of maximum in-plane Young’s modulus (fiber direction). Index 22 
denotes the direction transverse to the fiber in the plane of the lamina and index 33 denotes the direction 
perpendicular to the plane of the lamina. 


TABLE II. INPUT PARAMETERS. 


FE model 

global 

stabilization 

contact 

stabilization 

viscous 

regularization 

release 

tolerance 

last increment 

DCB-st 3 

2 10 -5 



0.2 

381 

DCB-st 4 

2 10 -6 



0.2 

1002 

DCB-st 5 

2 10 -7 



0.2 

550 

DCB-st6 

2 10' 8 



0.2 

1002 

DCB-st 7 

2 10 -8 



0.02 

1002 

DCB-st8 

2 10 -8 



0.002 

451 

DCB-ctl 


1 10 -5 


0.2 

1002 

DCB-ct 2 


1 10 -6 


0.2 

1002 

DCB-ct 3 


1 10 -7 


0.2 

751 

DCB-ct 4 


1 10 -7 


0.02 

1002 

DCB-ct 5 


1 10 -7 


0.002 

1002 

DCB-ct6 


1 10 -3 


0.002 

911 

DCB-vrl 



1 10 4 

0.5 

1002 

DCB-vr 2 



1 10 4 

0.3 

273 

DCB-vr 3 



1 10 -5 

0.5 

1002 

DCB-vr 4 



1 10' 5 

0.3 

1002 

SLB-stl 

2 10 -5 



0.2 

266 

SLB-st 2 

2 10 -6 



0.2 

1002 

SLB-ctl 


1 10 -6 


0.2 

133 

SLB-ct6 


1 10 -5 


0.5 

811 

SLB-ct8 


1 10 4 


0.5 

1002 

SLB-vrl 



1 10 4 

0.5 

65 

SLB-vr6 



1 10 -2 

0.5 

88 

SLB-vrl 2 



1 10 _1 

0.9 

537 
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Interlaminar sliding shear 
Mode II 

Figure 1. Fracture Modes. 



Interlaminar scissoring shear 
Mode III 
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0 experimental data 
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a. T300/914C. 
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mixed mode ratio G t /G T 
b. C12K/R6376. 

Figure 2. Mixed-mode fracture criterion for composite materials. 
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B 25.0 mm 

2 h 3.0 mm 

2 L 150.0 mm 

a 30.0 mm 


T300/914C [0] 24 


D±30: Cl 2K/R6376 [±30/0/-30/0/30/0 4 /30/0/-30/0/-30/30/^-30/30/0/30/0/ -30/0 4 /-30/0/30/±30] 


a. Double Cantilever Beam Specimen (DCB) 



B 25.4 mm 

t 1 2.032 mm 

t 2 2.032 mm 

2 L 177.8 mm 

a 34.29 mm 


D±30: C12K/R6376 [±30/0/-30/0/30/0 4 /30/0/-30/0/-30/30/^-30/30/0/30/0/ -30/0 4 /-30/0/30/±30] 

b. Single Leg Bending Specimen (SLB) 


Figure 3. Specimen configurations. 
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dimensions 



a. Deformed model ofDCB specimen with refined edges 



b. Deformed model of a DCB specimen for VCCT for ABAQUS analysis 
Figure 5. Full three-dimensional finite element models of a DCB specimen. 
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dimensions 




b. Deformed model of a DCB specimen with a coarse mesh 
Figure 6. Coarse full three-dimensional finite element models of a DCB specimen. 
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a. Deformed model of DCB specimen with a fine mesh 
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B 25.4 mm 

2 L 177.8 mm 

M 180.0 mm 

t 1 2.0 mm 
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a 34.0 mm 
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b. Deformed model of a DCB specimen with a coarse mesh 


dimensions 

B 25.4 mm 

2 L 177.8 mm 

M 180.0 mm 

t 1 2.0 mm 

t 2 2.0 mm 

a 34.0 mm 


z 




c. Deformed model of a DCB specimen with a coarse mesh 
Figure 7. Continuum-shell finite element models of a DCB specimen. 
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bonded 
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[ 0] 4 
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D±30: Cl 2K/R6376 [±30/0/-30/0/30/0 4 /30/0/-30/0/-30/30/^-30/30/0/30/0/ -30/0 4 /-30/0/30/±30] 
Figure 8. Deformed model ofDCB specimen and detail of region around delamination front. 
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the width of the specimen 
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detail of local 
model around 
delamination front 


B 



[30/0/-30/0] 
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30/30/0/30/0/ -30/0,/-30/0/30/±30] 


Figure 9. Deformed model ofSLB specimen and detail of region around delamination front. 
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D±30: Cl 2K/R6376 [±30/0/-30/0/30/0,/30/0/-30/0/-30/30/, -30/30/0/30/0/ -30/0,/-30/0/30/±30] 


Figure 10. Deformed model ofSLB specimen with refined center section. 
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Figure 11. Computed strain energy release rate distribution across the width of 
a DCB specimen ( model Figure 5a). 
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Figure 12. Failure index distribution across the width of a DCB specimen ( model Figure 5a). 
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load P, N 


-e — a=30mm 
-B — a=31mm 
-O — a=32mm 
-X — a=33mm 

—\ a=34mm 

-A — a=35mm 
-v — a=40mm 



0.4 0.6 0.8 1.0 

applied opening displacement (5/2, mm 


Figure 13. Load-displacement behavior ofDCB specimens with different delamination lengths a. 
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Figure 14. Benchmark: Critical load-displacement behavior for a DCB specimen. 
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load P, N 



Figure 15. VCCT for ABAQUS: Computed critical load-displacement behavior 
for DCB specimen obtained from results with global stabilization. 


load P, N 



applied opening displacement d/2, mm 

Figure 16. VCCT for ABAQUS: Computed critical load-displacement behavior 
for DCB specimen obtained from results with automatic global stabilization. 
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Figure 17. VCCT for ABAQUS: Computed critical load-displacement behavior 
for DCB specimen obtained from results with global stabilization. 




applied opening displacement 512, mm 

Figure 18. Computed critical load-displacement behavior for DCB specimen 
obtained from results with contact stabilization. 
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load P, N 



applied opening displacement 512 , mm 

Figure 19. VCCT for ABAQUS: Computed critical load-displacement behavior 
for DCB specimen obtained from results with contact stabilization. 
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Figure 20. Computed critical load-displacement behavior for DCB specimen 
obtained from results with viscous regularization. 
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delamination 
length 
a, mm 



Figure 21. Critical delamination length-displacement behavior for DCB specimen. 
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load P, N 

Figure 22. Computed critical delamination length-load behavior for DCB specimen. 
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Figure 23. Computed critical delamination length-displacement behavior for 
DCB specimen obtained from results with global stabilization. 
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Figure 24. Computed critical delamination length-load behavior for 
DCB specimen obtained from results with global stabilization. 
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load P, N 


energy, 

J 



applied opening displacement (5/2, mm 
Figure 25. Total energy calculation for DCB specimen. 



Figure 26. Calculated total strain energy in the model 
for DCB specimen ( model in Figure 7a). 
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Figure 27. Computed total strain energy and stabilization energy in the model of a DCB specimen 

obtained from results with global stabilization. 


energy, 

J 



Figure 28. Computed total strain energy and stabilization energy in the model of a DCB specimen 
obtained from results with automatic global stabilization. 
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Figure 29. Computed total strain energy and stabilization energy in the model of a DCB specimen 

obtained from results with global stabilization. 
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Figure 30. Computed total strain energy and stabilization energy in the model of a DCB specimen 

obtained from results with contact stabilization. 
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Figure 31. Computed total strain energy and stabilization energy in the model of a DCB specimen 

obtained from results with contact stabilization. 


energy, 

J 



Figure 32. Computed total strain energy and stabilization energy in the model of a DCB specimen 

obtained from results with viscous regularization. 
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initial delamination propagated fronts 

front propagation in 

teflon insert center first 

a. Scan of fractured DCB specimen 



b. Delamination front after 1000 increments ( Bottom surface ofFE model in Figure 4b) 
Figure 33. Delamination front shape for a DCB specimen. 
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Figure 34. Computed strain energy release rate distribution across the width of 

a multidirectional DCB specimen. 
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Figure 35. Failure index distribution across the width of a multidirectinal DCB specimen. 
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load P, N 



applied opening displacement <5/2, mm 

Figure 36. Load-displacement behavior of multidirectional DCB specimens 
with different delamination lengths a 



applied opening displacement S/2, mm 


Figure 37. Benchmark: Critical load-displacement behavior for a multidirectional DCB specimen. 
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applied opening displacement 612 , mm 

Figure 38. VCCT for ABAQUS: Computed critical load-displacement behavior 
for DCB specimen obtained from results with global stabilization. 
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Figure 39. Computed total strain energy and stabilization energy in the 
model of a multidirectional DCB specimen . 
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detail of local 


model around 
delamination front 


location of initial delamination front 


a. Delamination front obtained with contact stabilization ( Bottom surface of model in Figure 8) 



b. Detail of delamination front obtained with viscous regularization 
Figure 40. Delamination front shape for a multidirectional DCB specimen. 
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Figure 41. Computed strain energy release rate distribution 
across the width of a SLB specimen. 
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Figure 42. Distribution of mixed-mode ratio across the width of a SLB specimen. 
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Figure 43. Failure index distribution across the width of a SLB specimen. 
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Figure 44. Critical load-displacement behavior for SLB specimen. 
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Figure 45. Critical load-displacement behavior for SLB specimen. 



applied center deflection w , mm 

Figure 46. Computed critical load-displacement behavior for SLB specimen 
obtained from results with global stabilization. 
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Figure 47. Computed critical load-displacement behavior for SLB specimen 
obtained from results with automatic global stabilization. 



Figure 48. Computed critical load-displacement behavior for SLB specimen 
obtained from results with contact stabilization. 
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load P, N 



applied center deflection w, mm 

Figure 49. Computed critical load-displacement behavior for SLB specimen 
obtained from results with viscous regularization. 


load P, N 



applied center deflection w, mm 

Figure 50. Computed critical load-displacement behavior for SLB specimen 
obtained from results with stabilization. 
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Figure 51. Critical delamination length-displacement behavior for SLB specimen. 
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Figure 52. Computed critical delamination length-load behavior for SLB specimen. 
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Figure 53. Computed critical delamination length-displacement behavior for SLB specimen. 
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Figure 54. Computed critical delamination length-displacement behavior for SLB specimen. 
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Figure 55. Calculated total strain energy in the model of a SLB specimen ( model in Figure 9). 
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Figure 56. Computed total strain energy and stabilization energy in the model of a SLB specimen 

obtained from results with global stabilization. 
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Figure 57. Computed total strain energy and stabilization energy in the model of a SLB specimen 

obtained from results with automatic stabilization. 
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Figure 58. Computed total strain energy and stabilization energy in the model of a SLB specimen 

obtained from results with contact stabilization. 
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Figure 59. Computed total strain energy and stabilization energy in the model of a SLB specimen 

obtained from results with viscous regularization. 
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Figure 60. Computed total strain energy and stabilization energy in the model of a SLB specimen 

obtained from results with stabilization. 
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a. Delamination front obtained with contact stabilization (increment 76) 
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b. Detail of delamination front obtained with contact stabilization (increment 76) 
Figure 61. Computed delamination front shape for a SLB specimen ( model in Figure 9). 
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a. Delamination front obtained with contact stabilization ( increment 1000) 



b. Detail of delamination front obtained with contact stabilization (increment 1000) 
Figure 62. Computed delamination front shape for a SLB specimen ( model in Figure 9). 
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a. Delamination front obtained with global stabilization (increment 476) 



b. Detail of delamination front obtained with global stabilization ( increment 1000) 
Figure 63. Computed delamination front shape for a SLB specimen ( model in Figure 10). 
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